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Abstract 

Stress tensors are derived for the multiparticle collision dynamics algorithm, a particle-based 
mesoscale simulation method for fluctuating fluids, resembling those of atomistic or molecular sys- 
tems. Systems with periodic boundary conditions as well as fluids confined in a slit are considered. 
For every case, two equivalent expressions for the tensor are provided, the internal stress tensor, 
which involves all degrees of freedom of a system, and the external stress, which only includes the 
interactions with the confining surfaces. In addition, stress tensors for a system with embedded 
particles are determined. Based on the derived stress tensors, analytical expressions are calculated 
for the shear viscosity. Simulations illustrate the difference in fluctuations between the various de- 
rived expressions and yield very good agreement between the numerical results and the analytically 
derived expression for the viscosity. 
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I. INTRODUCTION 



Soft matter systems, such as colloidal suspensions or polymer and biopolymer solutions 
possess a wide range of length and time scales. The need to bride the length- and time-scale 
gaps for studies of these systems requires a simplified and coarse-grained description of the 
solvent degrees of freedom. Several mesoscale simulations techniques have been developed 
to meet this goal, which adequately reproduce fluid behavior. Among them, the multipar- 
ticle collision dynamics (MPC) method, originally proposed by Malevanets and Kapral,- 1 ^ 
has attracted considerable attention over the last few years. In a wide spectrum of applica- 
tions, it has been shown that MPC reproduces fluid properties adequately and accounts for 
hydrodynamic interactions, as illustrated in the recent review articles Refs. ^,4. 

Traditionally, there is a fundamental interest in the transport properties of complex flu- 
ids. The coarse-grained simulation approaches provide access to hydrodynamic phenomena 
on the mesoscale. It has been shown that MPC is very well suited to study non-equilibrium, 
rheological, and viscoelastic properties of such fluids.-^^*^ 1 ^ 1 ^' To fully characterize the 
equilibrium and non-equilibrium physical properties of the fluid system, adequate micro- 
scopic expressions — such as the stress tensor — have to be provided in order to establish a 
link between the simulation degrees of freedom and the macroscopic material properties, 
e.g., the viscosity. Particular expressions for the stress tensor of an MPC fluid have been 
provided in Refs. for a periodic system and in Ref. [l0| for a slit geometry. Analytical 

expressions for its viscosity have been derived by various approaches .^iMiIiiI£iI&iIZii& 

In this article, we will provide stress tensors at equilibrium and under shear flow for 
an MPC fluid as well as for a system with embedded point-like particles, which resembles 
the virial formulation of molecular systems.— i22i2Ii22i22i2ii2£i2& This formulation allows for a 
straightforward calculation of the stress and provides an expression for the solvent-solute 
coupling. Three dimensional systems with periodic boundary conditions are considered as 
well as fluids in a slit geometry, which requires an adaptation of the stress tensor due to 
wall interactions. Two equivalent formulations of the stress tensor are provided in every 
case,— >2ii2!i corresponding either to the mechanical definition of stress as force per area or 
as momentum flux across a hypothetical plane. — The instantaneous values of the respective 
expressions are different, but their averages are identical. Since the provided expressions 
are novel, we derive the shear viscosity from them, for both, a three dimensional periodic 



system as well as a system confined between walls under shear. We propose a modification 
of the MPC algorithm in the presence of walls with respect to the inclusion of wall-phantom 
particles. Compared to the original algorithm,— our formulation prevents any surface slip 
of fluid particles. 

The paper is organized as follows. In Sec. II, the multiparticle collision dynamics method 
is described as well as its coupling to a solute composed of mass points. In addition, the 
simulation parameters are listed. Stress tensors are determined for MPC fluids with periodic 
boundary conditions and those confined between two parallel walls without external field 
in Sec. III. The stress tensors in presence of shear flow for the same boundary conditions 
are calculated in Sec. IV. In Sec. V, analytical expressions for the viscosity are determined 
exploiting the derived stress tensors. Section VI summarizes our findings. Additional aspects 
of the fluid confined between surfaces are discussed in Appendix A, namely the center-of- 
mass velocity in a surface cell, and the wall collisional and wall kinetic stress tensors. 

II. THE MODEL 

A. Multiparticle Collision Dynamics 

In the MPC algorithm, a fluid is represented by point particles of mass m, which interact 
with each other by a stochastic process. The algorithm consists of alternating streaming 
and collision steps.-^ In the streaming step, the N s particle move ballistically and their 
positions change according to 

n(t) = ri(t-h)+hvi(t-h), (1) 

i — 1, . . . , N s , in the time interval h, which we denote as collision time. In the collision step, 
particles are sorted into cubic cells of side length a and their relative velocities with respect 
to the center-of-mass velocity of every cell are rotated around a randomly oriented axis by a 
fix angle a. This imposed stochastic process represents the effect of many real collisions. In 
a collision step, mass, momentum and energy are conserved which leads to the build up of 
correlations between the particles and gives rise to hydrodynamic interactions. Hence, the 
velocity of a particle changes according to 

Vi(t) = Vi(t) + (R(a) - E)(vi(t) - v cm (t)), (2) 



where Vi(t) is the velocity before the collision, R(a) is the rotation matrix,-^ v cm = 
Y2f c Vj = i/N c is the center-of-mass velocity of the particles contained in the cell of parti- 
cle i, and N c is the total number of fluid particles in that cell. E is the unit matrix. Hence, 
the change of momentum in a collision is 

Api(t) = m{vi{t) - Vi(t)) (3) 
= m(R(a)-E)[v i (t)-v an (t)]. 

Without external field, Vi(t + h) = Vi(t). In the presence of such a field, however, the 
velocity may change during the streaming step. Depending on the external field, additional 
forces have to be included in Eq. ([T]) .— »i2a2£ To insure Galilean invariance, a random shift 
is performed at any collision step.— Various alternative schemes for the stochastic process 
have been proposed by now.— ^ However, the actual collision process is not important for 
the derivation of a stress tensor, but it affects the dependence of the viscosity on the MPC 
parameters. 

B. Solute dynamics, solvent-solute coupling 

In complex fluids, solute particles are embedded in the MPC solvent. Here, we will assume 
that the solute is composed of mass points, e.g., polymers,— which interact with each other 
by pairwise potentials and their dynamics is treated by molecular dynamics simulations 
(MD). More complex objects, such as vesicles or solid bodies can also be embedded and their 
dynamics be coupled to the fluid.- We will consider N p objects (polymers) each composed 
of N m particles. The equations of motion of particle k with mass of object v read 

M k rl = F%, (4) 

where F£ is the total force. These equations are solved by, e.g., a velocity Verlet algorithm,— 
which provides the positions and velocities v%(t) starting at a time t — h. 

The solute particles can easily be coupled to the solvent by incorporating them in the 
collision step.— ^ For a collision cell with N c fluid particles and solute particles, which 



may belog to different objects, the center-of-mass velocity is given by 



Vein [tj 



i=l v k 



(5) 



N c m + J2 M k 



k=l 



which yields the momentum change (J2J) 



Apl{t) = M k {R(a) - E)[v%(t) - v cm (t)]. 



(6) 



Here, v and k belong to those polymers and monomers, respectively, which are within the 
considered collision cell. This results in an exchange of momentum between the solvent and 
solute degrees of freedom. The new monomer velocities are then used as initial conditions 
for the MD simulation of the embedded particles. Typically several MD steps are performed 
between multiparticle collisions, because the applied force fields require an integration time 
step, which is typically smaller than the collision time. 

C. Simulation parameters 

In a simulation, a cubic system is considered with linear extension L and an average 
number of N c = 10 particles in a collision cell. The rotation angle is set to a = 130°. 
Length and time are scaled according to fp = rp/a and t = t^JksT /ma 2 , which corresponds 
to the choice /c#T = 1, m = 1, and a = 1, where T is the temperature and ks the 
Boltzmann constant. The collision time h = 0.1 is applied, which is well in the collision 
dominated regime of the fluid dynamics.— 1 ^ In the calculation of the viscosity, shear is 
imposed either by Lees-Edwards boundary conditions^ or by the opposite movement of the 
confining parallel walls, with the shear rate 7 = 10~ 2 a/ k^T / ma? . 



III. STRESS TENSOR: NO EXTERNAL FIELD 



The actual form of the stress tensor depends on the boundary conditions of the fluid. 
Here, we will address periodic boundary conditions and solid walls. In general, the equation 
of motion of the ath (a G {x,y,z}) spatial component of the ith mass point is (r = 



TTliT'ia Fj a . (7) 

In case of periodic boundary conditions, T*j referrers to the position of the particle in the 
infinite system, i.e., we do not jump to an image, which is located in the primary box, when 
a particle crosses a boundary of the periodic lattice. Hence, 7*j is a continuous faction of 
time. Multiplication of Eq. (I7j) by and summation over all N s particles yields 

^ N a Ns Ns 

i=l i=l i=l 

The average over time (or an ensemble) yields 

N s \ / N s \ 

^2 m iVi a v i( i ) + ( ^ F ia r i/3 \ = 0, (9) 

because the term on the left hand side of Eq. (JHJ) vanishes for a diffusive or confined 
system.— 1 ^ Equation (Q will be the basis for the derivation of stress tensors. 

We will exploit the mechanical definition of the stress tensor given by <j a p = Fa/Ap, 
where F a denotes the total force in the spatial direction a across the surface of area Ap with 
normal in the spatial direction (3. 



A. Periodic boundary conditions 

For a system with periodic boundary conditions, we assume that initially all fluid particles 
are in the same box of the periodic system, which we will denote as primary box. In the 
course of time, the particles will diffuse out of that box. Some of them may reenter and leave 
again several times. The periodic images of particle i are located at the positions r\ + R n , 
with the lattice vectors 

Rn — ( n xL x ,n y Ly,n z L z ) T ', (10) 

corresponding to the lattice of images of the primary box. The n a s are integer numbers and 
L a denotes the box length along the a-direction. For a cubic lattice, L = L x = L y = L z = 



The potential energy of the solute particles comprises inter- and intramolecular pairwise 
contributions 

N p N p N m N m 

= 2 E E E E E u »m - < - 

v=l fi=l k=l 1=1 n 

N p N m N m 

+ 2EEE«- r n» ( n ) 

u=l k=l 1=1 

where the interaction of a particle with itself (not necessarily its image) is excluded and U v 
includes all intramolecular potentials of the object v. The sum over n accounts for all the 
images of a particular particle. The total force F% (j3J) of point k of object v is then given 
by 

N m 

^ = E^«-<) 

1=1 

Np N m 

+ EEE^TK-rf-i^). (12) 

H=l 1=1 n 

In general, infinite contributions of the intermolecular interactions have be taken into ac- 
count. For short-range interactions, however, only nearest images contribute significantly to 
the dynamics of a particle and we introduce a potential cut-off, which is chosen such that 
self-interactions of an object v are prevented. 22 



1. MP C fluid 

The stress tensor of the bare MPC solvent is obtained from Eq. Evidently, either 
fluid particles themselves or their images are in the primary box. Denoting the position 
(image or real) of a particle in the primary box by r[(t), the particle position itself is given 
by Vi(t) = r[{t) + Ri(t), where Ri(t) is the lattice vector at time t. The force exerted on 
the particle during the MPC collisions at times t q is 

oo 

F\(t) = ^2Api(t)S(t-t q ). (13) 

q=0 



The time average of Eq. (Q then reads 

1 f T 

{Fianp) = hm — / F ia (t)r ip (t) dt 

T^oo 1 J Q 
I N I ft q 

= }} m TrJ2h ^Pia(t)r^{t)S(t - t q ) dt 
1 N 



9=1 
N 



n^od Nh 

q=l 

' Ap ia (t q )r' i/3 (t q ) + Ap ia (t q )R i p(t q )) T , (14) 



where we introduced the average over collision steps 

\ 



1 - 



9=1 



We define now an instantaneous external stress tensor crfg 



</3 = VhJ2 ApiaRi ?- (16) 



Vh 

i=l 

Similarly, we introduce an instantaneous internal stress tensor by^-^ 

1 N s N s 
i=l i=l 

According to Eqs. (jUJ) and (fT4"l) . the averages of the two terms are equal, i.e., (<t%q)t — ( a ap)T- 
Hence, we obtain two equivalent expressions for the stress tensor. Equation (1161) corresponds 
to the mechanical definition as force per area (Rip ~ L) and Eq. (|17|) follows from the 
momentum flux across a surface. Correspondingly, the external stress tensor includes only 
force terms, i.e., collisional contributions, whereas the internal stress tensor comprises kinetic 
and collisional contributions. 

In general, the pressure follows from the stress tensor via the relation p = — cr aa /3. 

Figure [T] displays the dependencies of the averages (p l ) N , (p e ) N (cf. Eq. ( 1T5|) ) of the 
internal and external pressures on the number of collision steps, which yield the macroscopic 
pressure p = (p l ) T = (p e ) T in the limit iV — > oo. Evidently, both expressions approach the 
same limiting value for a large number of collision steps. The fluctuations of the average 
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FIG. 1: Internal {p l ) N (blue) and external (p e ) N (green, large fluctuations) pressure as a function 
of the number of collision steps N. The collision time is h = 0.1 and the time t = Nh. 

external pressure are larger, since the number of particles included in the pressure calculation 
are smaller as compared to the internal pressure. Moreover, the fluctuations of p e itself are 
larger and increase like the square root of t with time, because the fluid particles diffuse 
through the infinite periodic system.— The pressure is given by the kinetic contribution 
p = X^i=i 171 (^»)t /(^0 = NcksT. This follows from the fact that the momentum change 
in a collision cell is independent of the actual positions of the particles, hence Api and 7*j 
are uncorrelated and the collisional contributions to the internal pressure/stress vanish. 




2. MPC fluid and embedded particles 



For the case of embedded particles, Eq. (Q reads 

N 3 Np Nm 

= Y m ( V ^ V H3) T + Y Y Mk (vkJk^T 

1=1 u=l k = l 



i=l 

N p N n 



v=\ k=l 



v r v v 
kla[ r kf3 ~ r l/3l/ T 



!/=l k,l=l 
j N P N m 

+ O Y Y Y { F kTa[ r k, 



k/3 r l(3\/ T i 



v,^=l k,l=l n 



when we use Eqs. ( TT2|) and JT1, and t 



As described in Refs. 
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(18) 



ie abbreviation Fj£f = F$Jr k -r x - R n ). 
for the solute and with the same strategy as for the 



bare solvent, we obtain the following instantaneous stress tensors 

j Np Nm 

° e ^ = ^Y Y Y F kiX R np 

u,fi=l k,l=l n 



1 N s 1 N p N m 

+ VhY + VhYY A P** R h> ( 19 ) 



Vh^ v^/i 

1=1 f=l fe=l 

</3 = - y m *to«t/9 - y Y Y M ^kA 

1=1 !/=l fc=l 

i=l i/=l fc=l 



-J^ f J ' Tib 

juYY / '/ - ' 



r i> 

Ma L r fc/3 — ' l/3\ 

u=l k,l=l 



Nrr 

2V 

i 

-^EEE OT'fc - - RnA- (20) 

i^,At=l fc,/=l « 

Again, the averages are equal (cr?g) r = ( cr a /?)r- The solvent-solute coupling is captured in 
the terms with the momenta Ap v ka of the monomers as well as in those for the fluid momenta 
Ap ia of collision cells containing monomers. 



B. Confining walls 

We will now determine the stress tensors for an MPC fluid confined between two solid 
walls. The walls are parallel to the xy-pl&ne and periodic boundary conditions are applied 
along the x- and y-directions. The center of the reference system is located in the middle 
between the two walls, i.e., the wall positions are z w = ±L/2. The equations of motion 
of the fluid particles are then modified by the wall interactions. We will assume no-slip 
boundary conditions, which we realize by the bounce-back rule, i.e., the velocity of a fluid 
particle is reverted when it hits a wall (Vj — > —Vi)M 

The random shift perpendicular to the walls is implemented as follows. Without random 
shift, the most upper and lower border of the collision cells coincides with the respective wall. 
To enable a random shift, an additional layer of (empty) collision cells is added below the 
lower wall. In a random shift, the whole collision lattice is shifted in the positive ^-direction 
by a uniformly distributed displacement Az, with < Az < a, as illustrated in Fig. [2j The 
random shift typically leads to partially occupied cells at the walls, which in turn causes 
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FIG. 2: Illustration of the random shift and the distribution of particles in collision cells cut by 
the walls. The walls are located at z = ±L/2. Under shear, the walls move with the velocities 
u = ±7-L/2. The phantom particles are located in the centers of the truncated parts of the 
collision cell at (L + Az)/2 and — [L + a — Az)/2, respectively. They move with the velocities 
Up = j(L + Az)/2 and u p = —j(L + a — Az)/2, respectively. The dashed-dotted line indicates the 
linear velocity profile. 

a violation of the no-slip boundary condition under shear=. To restore no-slip boundary 
conditions, typically virtual particles are added to every cell cut by a wall and occupied by 
a number of particles N sc smaller than the average number of particles iV c , such that the 
average particle density is restored. However, this does not completely prevent slip, because 
the average center-of-mass position of all particles in a collision cell — including the phantom 
particle — does not coincide with the wall. In order to fully account for the no-slip boundary 
condition, we propose the following modification of the original approach. To treat a surface 
cell on the same basis as a cell in the bulk, i.e., the number of particles satisfies a Poisson 
distribution with the average iV c , we take fluctuations in the particle number into account by 
adding N sp particles to every cell cut by a wall such that (N sp + N sc ) = N c . The momentum 
P of a virtual particle is taken from the Maxwell-Boltzmann distribution with the variance 
er 2 = mN sp kBT and, at equilibrium, zero average. (The case of a shear flow is discussed in 
Sec. V B.) There are various ways to determine the number N sp . For a system with two 



parallel walls, we suggest to use the number of fluid particles in the surface cell cut by the 
opposite wall. The average of the two numbers is equal to iV c . Alternatively, N sp can be 
taken from a Poisson distribution with average N c accounting for the fact that there are 
already N sc particles in the cell. Collisions are then performed with all the particles in the 
cells. The center-of-mass velocity of the particles in a boundary cell is 

1 



' Ns, 



ma 



(21) 



m(N sc + N sp ) 

Naturally, this type of collisions will affect the external stress tensor. 

The total force on a fluid particle % comprises contributions from collisions among particles 
and collisions with the walls, i.e., 



q=0 q=0 

where t™ is the time at which the particle hits a wall and Apf = 
change. By averaging over a collision interval, the force term in Eq. ([9]) becomes 

i=i I h »=i 

T N s 



(22) 

2mv; its momentum 



i=l 



8 



z/3, 



with the Heaviside function 



G(x) 



1, x > 
0, x < 



(23) 



(24) 



By adding and subtracting the term E^ igbc Ap ia [0(r iz ) — ^]L5 z /3 in Eq. (123]) . we obtain the 

nal stress tensors 

— ^ APiaRipil ~ $z 



instantaneous external and internal stress tensors 
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(25) 



a 



aP 



N 3 I Ns 

1=1 1=1 



V 
L 

Vh 



6(r i2 ) - 



SzP, 



(26) 



idbc 



with r' iz = r iz , including both, the contributions by the surfaces as well as by the periodic 
boundary conditions. The surface contribution to the stress tensor in Eq. (1251) clearly shows 
that stress is force/area.— 1 ^ The last term of Eq. (1261) originates from the finite range of the 
fluid-surface interaction. In MPC, the non-locality of the fluid collisions is responsible for 
the extended range. The collision interaction of a particle with the surface is of zero range, 
whereas for a finite range potential, e.g., a Lennard- Jones potential, another term would 



appear as discussed in Ref. 



24 



Simulations yield a similar time dependence of the pressure as for the periodic system, 
which is display in Fig. [TJ 

We will not explicitly discuss the inclusion of a solute in the calculation of the stress tensor 
here. A detailed derivation of the stress tensors for mixed confined and periodic molecular 
systems is presented in Ref. |24j and the contribution of the solvent-solute interaction is 



identical to the terms presented in Eqs. (1X9]) and (1201) . 

IV. STRESS TENSOR: SHEAR FLOW 

The presence of shear flow alters some of the terms of the fluid stress tensors. Therefore, 
we will discuss this type of external field in more detail. 

In general, shear is applied in the x-direction and the gradient direction is along the 
z-axis. 



A. Periodic boundary conditions 

Again, we will discuss a system with periodic boundary conditions first. Because of the 
external field, the time average of the left hand side of Eq. ([8]) does not vanish anymore. 
Neglecting fluctuations for the moment, the velocity Vi X of the linear flow profile is v i X = ^ri Z) 
with the shear rate 7. The time average (d(vi X r iz )dt) is then given by lim^oo rf z (T)/T. 
Since a particle is diffusing along the gradient direction, rf z (T) ~ T and the average is 
finite. In order to arrive at a vanishing term, we subtract the derivative of the velocity 



profile d(jri Z )/dt = ^V{ Z from both sides of Eq. (J7J). This leads to the modified equation 

d Ns Ns 

— ^ m i v ix ~ ir iz )r iz = ^ m i v ix ~ ir iz )v iz 

i=l i=l 

N a N a 

+ ^2 FixTiz ~ i mVizTiz - ( 27 ) 

i=l i=l 

Evidently, the (time) average of the left hand side vanishes. Applying the definition of the 
time average (TH|) . the velocity terms on the right hand side read as 

{(v ix - in z )v iz ) = (v iz v' ix ) T + y (vl) T , 

(v iz r iz ) = - ((v iz + v iz )r iz ) T (28) 

in the stationary state. Note that Vi(t q ) is the velocity before the collision and Vi(t q ) that 
after the collision. Similar to the notation for the positions, v' ix denotes the velocity in the 
primary cell of the periodic system, i.e., v ix = v' ix + ■jRi X . (The particle velocities along 
the other spacial directions are identical for each periodic image.) The original expression 
({vix — ifi Z )vi Z ) T reduces to (v' ix Vi Z ) T) because the average (vi Z r' iz ) T vanishes. We like to 
point out that the change from Vi, 7*j to t>-, r[ corresponds to the application of Lees- 
Edwards periodic boundary conditions in non-equilibrium simulations of simple shear. For 
the sake of completeness, we emphasize that the time and ensemble average of the last term 
on the right hand side of Eq. (1271) is Yli^i m 7 ( v iz r iz) — m^NgD, where D is the diffusion 
coefficient of an MPC particle. 

We are now in the position to define instantaneous external and internal stress tensors as 

N s . N s 

° e xz = y^Yl ^VixRiz ~2Y^2 m ^ Viz + diz "> Riz > ( 29 ) 

i=l i=l 

1 Ns 'h Ns 1 Ns 

°xz = ~yYl - 7^7 Yl mV ^ z ~Vh^l A PixKz, (30) 

i=l i=l i=l 

which obey the relation (u xz ) T = (o'xz)t- The presence of the external field leads to addi- 
tional terms contribution to the stress tensors compared to the expressions (1161) and f|T7|) . 
The extra term in <j l xz results from the streaming dynamics and vanish in the limit h — > 0. 
Since a discrete time dynamics is fundamental for the MPC method, the collision time will 
always be finite. 

An example of the time dependence of the internal and external stress tensors, i.e., {cr xz } N , 
(<J XZ ) N , under shear is shown in Fig. [3j Both expressions approach the same limiting value 
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FIG. 3: Internal \<? XZ } N (blue) and external {<J% Z ) N (green, large fluctuations) stress tensor as 



function of the number of collision steps. The collision time is h = 0.1. At t = 0, the system is in 
a stationary state. 

for a large number of collision steps. The fluctuations of the external stress tensor component 
are again larger. 

B. Confining walls 

For the system described in Sec. Ill B, shear is imposed by the opposite movement of 
the confining walls with the velocities u = ^z w = ±jL/2. The structure of the external 
stress tensor (1251) is maintained. However, the wall momentum Ap™ x changes to Ap™ x = 
—2mVi Z + 2mu. In addition, the phantom particles of the partially filled surface cells possess 
a finite velocity. Hence, the momentum P x is now determined from the Maxwell-Boltzmann 
distribution of the same variance as before, but with the velocity 



As illustrated in Fig. [5J Az is the fraction of the cells truncated by the wall at z w = 
L/2. Correspondingly, a — Az is the fraction of the cells truncated by the opposite wall. 
In our description, the phantom particles are located in the centers (along the z-axis) of 
the truncated parts of a surface cells. The advantage of this approach over the previous 
implementation is that the center-of-mass velocity of a surface cell is equal to the velocity 
of the wall, as shown in Appendix A. This implies no-slip at the wall. 

For the internal stress tensor, the time average of the terms Vi X Vi Z over a collision interval 




(31) 



is modified. The integral now becomes 



2u 

v ix (t)v iz (t)dt = v ix (t q )v iz (t q ) - —v iz (t q )M q 



(32) 



for a particle which collides with a wall at t™ in the interval t q — h < t™ < t q and Ai*//i = 
1 — (t q — t^)/h. Evidently, the average over Vi Z is non-zero, because the relevant particles 
move always towards the respective surface. Hence, a xz becomes 



a, 



mv ix v iz + 
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(33) 



The other components of the stress tensor are obtained via Eq. fl26l) . 

Simulations confirm that the center of mass velocity of the particles interacting with 
phantom particles in the surface cells are indeed equal to the velocity of the respective 
surface. Moreover, the time dependent averages of the internal and external stress tensors 
(cr xz ) N , and {(7 XZ ) N are similar to those displayed in Fig. [31 



V. VISCOSITY 



The derived expressions for the stress tensors are independent of any particular collision 
rule. Transport coefficients such as the viscosity of a system, however, depend on the apply 
collision procedure. 

Analytical expressions for the viscosity of an MPC fluid have been derived by various 
approaches.-^ 1 ^^^^ Since the stress tensors of Eqs. f l29|) . (1301) . and (T33|) are novel, we 
will here derive the viscosity based on these expressions for the stochastic rotation version 
of MPC described in Sec. II. 

In simple shear flow with the velocity field v x = 'jz, the viscosity rj is related to the 
stress tensor via rj = o xz j^^ where the (macroscopic) stress tensor follows from o xz = 
( (7 L)t = { a xz)r- F° r an MPC fluid, the stress tensor is composed of a kinetic and collisional 
contribution,-^^^^ i.e, a xz = a^ z + a x ° z , which implies that the viscosity r] = 77km + f]co\ 
consists of a kinetic 77km and collisional 77 co i part too.-^^ 1 ^^ 



A. Periodic boundary conditions 



For a system with periodic boundary conditions, the two contributions to the viscosity 
are conveniently obtained from the internal stress tensor (13"Uj) . 

The kinetic contribution ?7k in is determined by the streaming step, i.e., velocity dependent 
terms in Eq. (|30|) . To find the mean (v' ix v' iz ), we consider a complete MPC dynamics step. 
The velocity v' ix (t q ) before streaming is related to the velocity v' ix {t q + h) after streaming via 
v'ixttq + h) = v ix (t q + h) - jr iz (t q + h) = v ix (t q ) - jr iz (t q ) - jv iz (t q )h = v' ix (t q ) - jv iz (t q )h. 
With v' iz (t q + h) = v iz (t q ), we obtain the average 

(v' lx (t q + h)v iz (t q + h)) = (v' lx (t q )v lz (t q )} - ih (vl) . (34) 

Here, the average comprises both, a time average and an ensemble average over the orien- 
tation of the rotation axis. The velocities after streaming are changed by the subsequent 
collisions, which yields, with the corresponding momenta of the rotation operator H(a), 
(v' ix (t)v iz (t)) = f (v' ix (t)v iz (t)) and/ = l + (l-l/AT c )(2cos(2a) + 2cos(«)-4)/5.i^Note, 
velocity correlations between different particles are neglected, i.e., molecular chaos is as- 
sumed. Thus, in the steady stead [(v[ x {t)v' iz {t)) = (v' ix (t + h)v' iz {t + h))], we find 

(v' ix v iz ) N = - fzrj ( u te> ( 35 ) 

by using Eq. fl34l) . Hence, with the equipartition of energy (vf z ) = k B T/m, the kinetic 
viscosity is given by 
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V L(iV c -l)(4-2cos(a) -2cos(2a)) 2 
in agreement with previous calculations. 

The collisional viscosity is determine by the momentum change of the particles during 
the collision step. Since the collisions in the various cells are independent, it is sufficient to 
consider one cell only. The positions of the particles of that cell can be expressed as r\ = 
r c + Ari, where r c is chosen as the center of the cell. Because of momentum conservation, 
the term Y^/iti ^Vix r \ z then reads as Yli=i ^PixAr iz . The averages over thermal fluctuations 
and random orientations of the rotation axis yield 

(Ap ix Ar iz ) = ^p-(cos(a) - 1) 
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(37) 



The average over the uniform distribution of the positions within an cell yields (Ar^Ar^) = 
for i j and 

l r a/2 a 2 

- / Arid* = Yo- ^ 

a J -a/2 11 

Hence, the collisional viscosity is given by 

N s ma 2 , ., / 1 \ 

'^iw' 1 -^'! 1 "^ (39) 

again in agreement with previous calculations. 

Here, we assume that the number of particles in a collision cell N c is sufficiently large 
(N c > 3) to neglect fluctuations.- For a small number of particles, density fluctuations have 
to be taken into account. Then, Eqs. ( 1361) and (159"]) have to be averaged over the particle 
number using a Poisson distribution with the mean value N s /V.—^ 

We perform simulations for various MPC parameters and found a very good agreement 
between the viscosities determined via Eqs. ([29]) . fl30l) and the analytical expression Eqs. 
and (1391) . 



B. Confining walls 

Under confinement, the component a e xz of the external stress tensor is determined by 
the collisions of the fluid particles with the walls, which corresponds to the kinetic contri- 
bution, and the collisions of fluid particles within the partially filled surface cells, which 
yields the collisional contribution to the viscosity. Since the averages over the stress tensor 
contributions from each wall are equal, we find 

I „ e \ ~kin , „col 
(°W =°xz + a xz 

\i=l / \i£bc I 

and the averages are taken over one surface only. 

As shown in Appendix B and C, the evaluation of the averages yields exactly the same 
expressions for the viscosities as derived for a periodic system in Sec. IV A, namely Eqs. 
( 1361) and ( |39l) . We like to point out that this is not true in general. A simulation study with 
the mean of the momentum P x = m(N c — N sc )u yields different results for the two viscosity 
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FIG. 4: Viscosities determined via the internal (bullets) and external (open squares) stress tensors 
for a system confined between walls as function of the collision time. The analytical results for the 
total (black), the kinetic (red, ~ h), and collisional (blue, ~ l/h) contributions are presented by 
solid lines. 

contributions, although the total viscosity agrees with the theoretical prediction. 10 Only for 
our choice of the mean momentum (13T1) follows agreement with the theoretical expressions. 

Figure 0] depicts viscosities determined via the internal {cr l xz ) T fl33|) and external (cr^ 2 ) T 
( 140]) stress tensors and their respective collisional and kinetic contributions. Evidently, the 
averages agree very well with each other. Moreover, the simulation results agree very well 
with the analytical predictions for the kinetic and collisional viscosities of Eqs. ( 1361) and 

dSHD- 

VI. SUMMARY 

In this article, we have been introducing external and internal stress tensors for an MPC 
fluid resembling those of atomistic molecular fluids. Systems with periodic boundary con- 
ditions and fluids confined in a slit have been addressed and their peculiarities have been 
worked out. Moreover, the modifications of the stress tensors caused by the presence of 
simple shear have been determined. Based on the derived stress tensors, an analytical ex- 
pressions for the viscosity has been derived, which agrees with previous results.-^ 1 ^^^ 
In addition, stress tensors for systems containing solute molecules are presented, which are 
coupled to the solvent in the MPC collision step. These expressions explicitly comprise the 
solvent-solute contributions to the stress tensors. 



The stress tensors can easily be modified to account for a different coupling between 



the solvent and the solute. In Refs. 



371 . the solute interacts through an intermolecular 



potential, i.e., the Lennard- Jones potential, with the solvent. This results in an additional 
virial term in the stress tensor — similar to the solute intermolecular interactions — with the 
forces between the solvent and the solute particles and their respective positions. 

Simulations for various MPC parameters confirm the equivalence of time averages of the 
internal and external stress tensors of the fluid for both types of boundary conditions. More- 
over, the calculated viscosities are in accord with the corresponding analytical expressions. 

The stress tensors can easily be calculated, since they require known quantities, i.e, 
positions, velocities, and momenta changes, only. Moreover, all particles contribute in 
the calculation of the internal stress tensors and no extra hypothetical plane needs to be 
introduced.— 1 ^ 
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APPENDIX A: CENTER-OF-MASS VELOCITY IN SURFACE CELLS 



The center-of-mass velocity of all particles in a surface cell truncated by a wall and filled 
with a phantom particle of momentum P is given by Eq. (l2~Tj) . The average of the component 
in the flow direction for the wall at z w = L/2 reads 
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where N c = N sc + N sp . The fluctuations of P x have been averaged out already. The average 
over the fluctuations of the fluid particle velocities yields 
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a — Az is the part of the intersected surface cell which is within the fluid slit, and the average 
z of the particle position in a cell is 

1 f L/2 L a-Az . . 

z = — / zdz = — . (A3) 

a — Az J L /2-a+Az 2 2 

The remaining average is over the random shift Az and the particle number N sc . The 
average of Az, < Az < a, yields (Az) = a/2 and (N sc ) = N c /2. Thus, we find (v cm ^ x ) = u. 



APPENDIX B: SURFACE COLLISIONAL STRESS TENSOR 

The collisional contribution to the stress tensor (HOI) at a wall is given by 



1 / Nsc \ 



= (Bi) 



because the contributions from the various cells are independent. Averaging over the orien- 
tation of the rotation axis and the fluctuations of the momenta yields 
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with Eq. (1A3I) . The number of particles N sc and N sp are binomially distributed,— which 
yields the average (N sp N sc ) = N C (N C — l)Az(l — Az/a)/a. Az/a and 1 — Az/a are the 
probabilities to find a particle in one of the respective parts of a collision cell. The average 
over the random shift Az (0 < Az < a) yields 

A2 fi_^\ = I rUl-*) dA z =°. (B3) 
V a)l a J V a I 6 



Thus, 



a Tz = Tpli 1 ~ c °s(«)) (A c - 1) (B4) 
loan 



and the collisional viscosity is given by Eq. (1391) 



APPENDIX C: SURFACE KINETIC STRESS TENSOR 



The kinetic contribution to the stress tensor (140]) at a wall is given by 

The average contains the information about the number of particles colliding with a wall 
in a streaming step, which can be determined by applying kinetic theory The number of 
particles in a volume element Vdzdv/L of the one-particle phase-space is given by dN = 
N s P(y)dzdv/L, with P(v) the velocity distribution function. Hence, Eq. (1C1I) can be 
reformulated as 

AT r j-L/2 

= -£ / / A P :P(v) dzdv (C2) 

11 J JL/2-hv z 

for the surface at z w = L/2. Only particles with velocities v z h > are able to reach 

the surface in the collision time interval h. By substitution of the momentum and with 

v' x = v x — jz, we obtain 

kin 2mN c f f L/2 . x J J 

^ = - —r^ / / fa- u)P(v) dzdv 

n J JL/2-hv z 

2mN r f f L/2 . , . N „, , 7 , 

(v x — it + 7^J-r (t>J a^ai; 

*L/2~hv z 



h 
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Extending the velocity integration from < w z < oo to — oo < v z < oo yields 

„kin 
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The velocities in this equation are the velocities after collision. With Eq. (1341) . the above sta- 
tionary state correlation function can be replaced by the correlation function of the velocities 
after streaming, which gives 

„kin 



-mN c ( (v'J z ) + ^ {vl) ) . (C5) 



2 

This equation agrees with the corresponding expression in Eq. (15U1) and we therefore obtain 
the same analytical expression for the kinetic viscosity as for the periodic system, namely 
Eq. (J36D. 
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